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We  have  analyzed  the  transferability  of  a  previously  proposed  intermolecular 
potential  for  nitramine  crystals  to  reproduce  the  experimentally  determined  crystal 
structures  (within  the  approximation  of  rigid  molecules)  of  51  nitro  compounds.  These 
compounds  include  different  types  of  acyclic,  monocyclic,  and  polycyclic  molecules.  It  is 
shown  that  this  potential  model  accurately  reproduces  the  experimentally  determine 
crystallographic  structures  and  lattice  energies  for  the  majority  of  these  crystals.  Further 
testing  of  the  proposed  intermolecular  potential  has  been  done  by  performing  isothermal- 
isobaric  molecular  dynamics  (MD)  simulations  over  the  temperature  range  100-450  K,  at 
atmospheric  pressure,  for  the  monoclinic  phase  of  the  2,4,6-trinitrotoluene  (TNT)  crystal 
and  for  the  polymorphic  phase  I  of  the  pentaerythritol  tetramtrate  (PETN  I)  crystal.  In 
each  case,  the  results  show  that  throughout  the  MD  simulations,  the  average  structures  of 
the  crystals  maintain  the  same  space  group  symmetry  as  the  one  determined 
experimentally,  and  there  is  a  good  agreement  between  the  calculated  crystallographic 
parameters  and  the  experimental  values. 
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We  have  analyzed  the  transferability  of  a  previously  proposed  intermolecular  potential  for  nitramine  crystals 
to  reproduce  the  experimentally  determined  crystal  structures  (within  the  approximation  of  rigid  molecules) 
of  51  nitro  compounds.  These  compounds  include  different  types  of  acyclic,  monocyclic,  and  polycyclic 
molecules.  It  is  shown  that  this  potential  model  accurately  reproduces  the  experimentally  determined 
crystallographic  structures  and  lattice  energies  for  the  majority  of  these  crystals.  The  best  agreement  with 
experimental  structural  and  energetic  data  is  obtained  when  the  electrostatic  charges  have  been  determined 
using  ab  initio  methods  that  include  electron  correlation  effects,  namely  MP2  and  B3LYP.  The  use  of  the 
electrostatic  charges  calculated  at  the  Hartree— Fock  level  results  in  large  differences  between  the  predicted 
and  the  experimental  values  of  the  lattice  energies.  This  difference  can  be  significantly  decreased  by  scaling 
the  electrostatic  charges  with  a  general  factor  without  introducing  significant  variations  of  the  predicted 
crystallographic  parameters.  Further  testing  of  the  proposed  intermolecular  potential  has  been  done  by 
performing  isothermal-isobaric  molecular  dynamics  (MD)  simulations  over  the  temperature  range  100-450 
K,  at  atmospheric  pressure,  for  the  monoclinic  phase  of  the  2,4,6-trinitrotoluene  (TNT)  crystal  and  for  the 
polymorphic  phase  I  of  the  pentaerythritol  tetranitrate  (PETN  I)  crystal.  In  each  case,  the  results  show  that 
throughout  the  MD  simulations  the  average  structures  of  the  crystals  maintain  the  same  space  group  symmetry 
as  the  one  determined  experimentally  and  there  is  a  good  agreement  between  the  calculated  crystallographic 
parameters  and  the  experimental  values.  The  thermal  expansion  coefficients  calculated  using  the  present  model 
indicate  an  overall  anisotropic  behavior  for  both  TNT  and  PETN  I,  with  a  thermal  isotropy  for  PETN  I  along 
cell  directions  a  and  b. 


I.  Introduction 

The  work  presented  here  is  the  fifth  in  a  series  of  studies1-4 
that  investigates  the  degree  to  which  an  intermolecular  potential 
energy  function  that  was  developed  to  describe  a  single 
molecular  crystal  can  be  extended  to  describe  crystal  structures 
of  other  similar  systems.  The  function,  consisting  of  atom- 
atom  (6-exp)  Buckingham  terms  and  electrostatic  interaction 
terms  in  the  form  of  partial  charges  associated  with  the  atoms, 
was  parametrized  to  reproduce  the  experimental  crystal  structural 
information  of  the  a-form  of  the  solid  explosive,  RDX  (hexahy- 
dro-l,3,5-trinitro-l,3,5-.s-triazine).1  The  parametrization  of  the 
function  was  done  such  that  molecular  packing  calculations 
(MP)  reproduced  the  experimental  structure  of  the  crystal  and 
its  lattice  energy  with  the  electrostatic  charges  determined  at 
the  second-order  Moller— Plesset  6-3 1G**  level.  This  intermo- 
lecular  potential  was  also  used  in  isothermal— isobaric  molecular 
dynamics  (NPT-MD)  calculations  at  ambient  pressure  for 
temperatures  ranging  from  4.2  to  325  K.  The  results  of  the 
simulations  are  in  good  agreement  with  experiment,  with  the 
lattice  dimensions  being  within  2%  of  experiment  and  almost 
no  rotational  or  translational  disorder  of  the  molecules  in  the 
unit  cell.  The  space  group  symmetry  was  maintained  throughout 
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the  simulations.  Thermal  expansion  coefficients  were  also 
determined  for  the  model  and  are  in  reasonable  accord  with 
experiment. 

The  recently  developed  explosive  2,4,6,8, 10,12-hexanitro- 
hexaazaisowurtritane  (HNIW)  can  be  described  as  a  bridged 
pair  of  RDX  molecules,  suggesting  that  the  intermolecular  forces 
for  HNIW  might  be  similar  to  those  of  RDX.  To  explore  this 
possibility,  we  performed  MP  and  NPT-MD  simulations2  for 
three  of  the  polymorphs  of  HNIW  (/?-,  e-,  and  y-phases)  at 
ambient  pressure  and  over  the  temperature  range  4.2  to  425  K 
using  the  form  of  potential  used  in  the  RDX  study.1  The 
parameters  for  the  Buckingham  terms  remained  unchanged,  and 
the  Coulombic  interaction  terms  between  electrostatic  charges 
were  determined  from  fits  to  ab  initio  electrostatic  potentials 
calculated  for  the  individual  molecules  corresponding  to  the 
different  polymorphs,  whose  atoms  are  arranged  in  the  experi¬ 
mental  configurations.  We  found  that  the  potential  predicts  the 
right  order  of  stability  for  different  phases  of  HNIW  (c  >  > 

y)  crystals  in  agreement  with  experimental  measurements.4  At 
300  K,  the  average  lattice  dimensions  agree  very  well  with 
experimental  values,  with  the  corresponding  differences  for  the 
individual  cell  edge  lengths  being  no  more  than  1.0%  for  the 
e-polymorph,  0.9%  for  the  ^-polymorph,  and  2.5%  for  the 
y -polymorph.  For  the  e-  and  y-phases,  the  variations  of  the  unit 
cell  angle  from  the  experimental  values  are  1.3%  and  0.1%, 
respectively,  while  the  other  two  angles  of  the  unit  cell  remain 
approximately  equal  to  90°.  For  the  /?~phase,  all  three  crystal- 
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lographic  angles  remain  approximately  equal  to  90°,  in  agree¬ 
ment  with  experiment.  Additionally,  little  rotational  or  trans¬ 
lational  disorder  occurs  throughout  the  simulations.  The  largest 
deviation  between  experimental  values  and  predictions  of 
molecular  orientation  occurs  for  y-HNIW;  the  predicted  value 
of  one  of  the  Euler  angles  defining  the  molecular  orientation 
deviates  by  4.4°  from  the  experimental  value. 

We  next  performed  MP  and  NPT-MD  simulations  at  atmos¬ 
pheric  pressure  and  different  temperatures  for  the  fi-,  a-,  and 
d-phases  of  another  nitramine  crystal,  the  explosive  1 ,3,5,7- 
tetranitro- 1 ,3,5,7-tetraazacyclooctane  (HMX).3  Again,  die  Buck¬ 
ingham  terms  were  those  used  in  the  RDX  study,  and  the  only 
differences  in  the  potential  are  the  Coulombic  interactions 
between  electrostatic  charges  centered  on  the  atoms.  At  room 
temperature,  the  predicted  average  lattice  dimensions  of  the 
£-phase  are  within  0.7%  of  the  experimental  values,  and  at  376 
K,  the  differences  of  the  lattice  dimensions  for  the  a-phase  are 
within  2.6%.  The  differences  for  the  (3-phase  are  within  4.4% 
at  433  K.  In  addition,  for  all  three  phases,  the  angles  of  the  unit 
cell  remain  close  to  the  experimental  values  with  the  maximum 
difference  of  2.5%  for  the  /?  angle  of  /?-HMX.  However,  it 
should  be  pointed  out  that  in  the  cases  of  a-  and  <5-HMX,  the 
crystallographic  data  were  determined  at  T  =  295  K,  although 
this  temperature  is  outside  the  stability  range  of  these  two 
phases.  There  are  no  significant  displacements  of  the  molecular 
center  of  mass  or  increase  of  the  degree  of  rotational  disorder 
for  the  three  phases.  The  largest  difference  between  experiment 
and  predictions  of  molecular  orientation  is  for  one  of  the  Euler 
angles  for  a-HMX;  the  predicted  average  value  is  3.8°  larger 
than  the  experimental  value.  Besides  the  geometrical  parameters, 
the  calculated  lattice  energies  for  the  three  phases  support  the 
experimentally  determined  polymorph  stability  ranking  (Ji  >  a 
>  (5)  given  by  McCrone.6  Moreover,  for  the  £-  and  d-phases, 
where  experimental  values  for  the  heats  of  sublimations  have 
been  determined,  the  predicted  lattice  energies  are  in  very  good 
agreement  with  the  experimental  values. 

More  recently,  we  have  extended  our  investigations  of  the 
transferability  properties  of  this  interaction  potential  to  a 
collection  of  30  nitramine  crystals.4  The  crystals  are  composed 
of  monocyclic,  polycyclic,  and  acyclic  nitramine  molecules.  The 
molecules  associated  with  the  nitramine  crystals  were  chosen 
as  representative  examples  of  acyclic  and  cyclic  nitramines.  In 
the  latter  case,  we  have  included  different  types  of  mono-  and 
polycyclic  nitramines,  particularly  crystals  of  importance  in 
energetic  materials.  For  most  of  the  crystals,  the  predicted 
structural  lattice  parameters  differ  by  less  than  2%  from  the 
experimental  structures  with  small  rotations  and  practically  no 
translations  of  the  molecules  in  the  asymmetric  unit  cell. 

In  the  present  study,  we  extend  our  investigations  beyond 
the  case  of  nitramine  crystals.  For  this  purpose,  we  have 
performed  MP  calculations  on  51  crystals  comprising  a  wide 
variety  of  compounds  such  as  nitroalkanes,  nitroaromatics, 
nitrocubanes,  polynitroadamantanes,  polynitropolycyclounde- 
canes,  polynitropolycyclododecanes,  hydroxynitro  derivatives, 
nitrobenzonitriles,  nitrobenzotriazoles,  and  nitrate  esters  such 
that  a  comprehensive  test  to  our  potential  can  be  achieved.  We 
have  been  particularly  interested  to  see  if  the  geometrical  and 
the  known  energetic  parameters  for  these  types  of  crystals  can 
be  reproduced  accurately  by  the  proposed  intermolecular 
potential. 

As  in  the  preceding  studies,1"4  we  used  the  RDX  Buckingham 
potential  plus  Coulombic  interactions  terms  obtained  through 
fitting  of  partial  charges  centered  on  each  atom  in  the  experi¬ 
mental  arrangement  of  the  molecule  to  a  quantum  mechanically 


derived  electrostatic  potential.7  Moreover,  as  in  the  case  of  the 
nitramine  crystals  4  we  have  investigated  how  the  geometrical 
and  energetic  parameters  predicted  in  molecular  packing 
calculations  depend  on  charges  calculated  from  ab  initio  methods 
that  do  or  do  not  include  electron  correlation  effects.  Specifi¬ 
cally,  we  used  different  sets  of  charges  derived  from  the 
Hartree— Fock  (HF)  wave  function8  or  from  methods  that 
employ  electron  correlations  such  as  second-order  Moller- 
Plesset  (MP2)9  and  density  functional  theory  (DFT).10  We  again 
note  that  the  main  limitation  of  the  calculations  is  the  assumption 
of  rigid  molecules,  but  this  model  can  be  used  to  study  processes 
at  temperatures  and  pressures  where  molecular  deformations 
are  negligible.  Our  intent  is  to  extend  the  model  to  allow  for 
deformation  of  the  molecules  by  incorporating  intramolecular 
interaction  terms. 

The  organization  of  the  paper  is  as  follows.  In  section  II,  the 
intermolecular  potential  used  to  simulate  the  crystals  is  pre¬ 
sented.  In  sections  HI,  the  details  of  calculations  using  molecular 
packing  methods  and  isothermal  -isobaric  MD  calculations  are 
described.  The  results  of  these  calculations  are  given  in  section 
IV.  The  main  conclusions  are  summarized  in  section  V. 

n.  Intermolecular  Potential 

In  this  work,  we  adopt  the  same  general  principles  for  atom- 
atom  potentials  that  proved  to  be  successful  in  modeling  of  the 
nitramine  crystals.1"4  In  particular,  we  assume  that  (a)  the 
intermolecular  interactions  depend  only  on  the  interatomic 
distances,  (b)  the  interaction  potential  can  be  separated  into 
contributions  identified  as  van  der  Waals  and  electrostatic,  and 
(c)  the  same  type  of  van  der  Waals  potential  is  used  for  the 
same  type  of  atoms,  independent  of  their  state  of  valence. 
Moreover,  we  assume  the  transferability  of  the  potential 
parameters  between  similar  molecules;  that  is,  we  extend  the 
validity  of  the  potential  parameters  determined  for  RDX  crystal 
to  all  the  nitro  compounds  considered  in  the  present  database. 

We  approximate  the  intermolecular  interactions  between  the 
molecules  of  the  crystal  as  a  pairwise  sum  of  Buckingham  (6- 
exp)  (repulsion  and  dispersion)  and  Coulombic  (C)  potentials 
of  the  form 

exp (-B^r)  -  C^r6  (1) 

and 


where  r  is  the  interatomic  distance  between  atoms  a  and  /?,  qa 
and  qp  are  the  electrostatic  charges  on  the  atoms,  and  e0  is  the 
dielectric  permittivity  constant  of  vacuum. 

The  parameters  for  the  6-exp  potential  in  eq  1  are  those 
previously  determined  for  the  RDX  crystal.1  We  use  the  same 
combination  rules  for  calculating  the  heteroatom  parameters 
from  homoatom  parameters  as  previously  reported.1  The  as¬ 
signments  of  the  electrostatic  charges  were  made  by  fitting  a 
*  set  of  atom-centered  monopole  charges  for  the  isolated  molecule 
to  reproduce  the  quantum  mechanically  derived  electrostatic 
potential,  which  is  calculated  over  grid  points  surrounding  the 
van  der  Waals  surface  of  the  molecules.  This  method  of  fitting 
the  electrostatic  potential  was  proposed  by  Breneman  and 
Wiberg7  and  is  incoiporated  in  the  Gaussian  94  package  of 
programs11  under  the  keyword  CHELPG  (electrostatic-potential- 
derived  atomic  charges).  The  quantum  mechanical  calculations 
have  been  done  at  the  Hartree— Fock  (HF),8  second-order 
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Mdller— Plesset  (MP2),9  and  density  functional  theory  (DFT)10 
levels  to  investigate  the  effect  of  electron  correlation.  The 
density  functional  that  was  used  includes  the  exchange  func¬ 
tional  described  by  the  fitted  three-parameter  hybrid  of  Becke12 
and  the  correlation  functional  of  Lee,  Yang,  and  Parr  (Becke3- 
LYP).13  All  of  the  above  theoretical  calculations  were  done  using 
the  basis  set  6-3 1G**  (split- valence  plus  d-type  and  p-type 
polarization  functions).14 

It  has  been  previously  shown1516  that  the  neglect  of  electron 
correlation  in  self-consistent  wave  functions  results  in  an 
overestimation  of  the  electrostatic  interactions  and  that  this 
overestimate  is  mainly  a  scaling  effect  A  scaling  factor  of  0.9 
was  shown  to  give  improved  agreement  between  the  calculated 
and  the  experimental  dipole  moments  for  a  set  of  eight  small 
molecules,15  in  a  study  of  the  electrostatic  interactions  of  a 
dipeptide,16  and  in  determining  the  crystal  structures  of  polar 
organic  molecules.17  Our  previous  study  of  the  nitramine 
crystals4  showed  that  the  use  of  the  0.9  scaling  factor  for  the 
electrostatic  charges  determined  at  the  HF  level  significantly 
improves  the  accuracy  of  the  predicted  lattice  energies  of  the 
crystals.  We  further  investigate,  in  the  present  study  different 
electrostatic  models  and  evaluate  the  effects  of  the  scaling 
procedure.  Specifically,  four  electrostatic  models  were  tested 
for  each  of  the  51  crystals.  Two  of  them  use  electron  correlation 
methods,  namely  MP2  and  B3LYP,  the  third  one  uses  unsealed 
HF  charges,  and  the  last  employs  the  HF  charges  scaled  by  0.9 
(denoted  as  0.9HF). 

HL  Computational  Approach 

Molecular  Packing  Calculations.  Molecular  packing  cal¬ 
culations  are  used  to  test  empirical  or  semiempirical  intermo¬ 
lecular  potential  energy  functions  of  organic  crystals.18’19  The 
calculations  minimize  the  lattice  energy  of  the  models  with 
respect  to  the  structural  degrees  of  freedom  in  the  crystals.  For 
crystals  in  which  the  asymmetric  unit  contains  one  molecule 
that  occupies  an  arbitrary  position,  the  maximum  number  of 
degrees  of  freedom  is  12.  These  correspond  to  the  six  unit  cell 
constants  (ay  b,  c,  a,  /?,  y),  the  three  rotations  ( 0\ ,  #2,  and 

the  three  translations  (n,  T2,  r3)  of  the  rigid  molecule.  The 
number  of  structural  degrees  of  freedom  might  be  reduced 
depending  on  the  symmetry  restrictions  of  the  space  group. 
Crystals  in  which  the  asymmetric  unit  contains  more  than  one 
molecule  have  additional  degrees  of  freedom  to  describe  the 
rotation  and  translation  of  the  molecules.  As  in  the  case  of 
nitramine  crystals 4  we  consider  that  the  crystals  can  be 
represented  as  an  ensemble  of  rigid  molecules. 

A  stable  crystal  configuration  is  obtained  by  assuming  that 
the  crystal  energy  can  be  represented  as  a  function  of  the 
structural  lattice  parameters  and  minimizing  the  crystal  energy 
with  respect  to  the  structural  lattice  parameters.  The  minimiza¬ 
tion  is  performed  using  traditional  steepest-descent  and  New¬ 
ton— Raphson  procedures.20*21 

Two  series  of  MP  calculations  were  performed.  In  the  first 
series,  the  energy  minimizations  were  performed  for  all  the 
crystals  using  the  program  PCK91.22  Starting  configurations 
correspond  to  the  experimentally  observed  geometries.  This 
program  employs  the  accelerated  convergence  method1*20  for 
accurate  evaluation  of  the  crystal  Coulombic  and  dispersion 
lattice  sums,  with  the  first  and  second  derivatives  of  the  crystal 
energy  evaluated  analytically.  The  space  group  symmetry  is 
maintained  throughout  the  energy  minimization,  reducing  the 
number  of  independent  variables  in  the  minimization  procedure. 
For  example,  nitromethane  (see  entry  1  in  Table  IS  of  the 
Supporting  Information)  has  space  group  symmetry  P2i2i2i  (Z 


=  4).  In  this  case,  the  crystallographic  parameters  varied  in  the 
minimization  using  the  PCK91  program  are  the  three  dimensions 
of  the  unit  cell  and  the  three  rotations  and  translations  of  the 
molecule  in  the  asymmetric  unit  cell.  The  three  angles  a,  /?, 
and  y  of  the  unit  cell  were  set  at  90°  and  were  not  allowed  to 
vary.  The  structural  shift  factor  2h23 

F  =  (A0/2)2  +  (IOAjc)2  +  (100Aa/a)2  +  (100AM>)2  + 

(100Ac/c)2  +  (Aa)2  +  (Apf  +  (Ay)2  (3) 

provides  a  measure  of  the  quality  of  the  predicted  geometrical 
crystallographic  parameters  relative  to  the  experimental  values; 
A 6  is  the  total  root-mean-square  (mis)  rigid-body  rotational 
displacement  (in  degrees)  after  minimization.  Ax  is  the  rms  total 
rigid-body  translational  displacement  (in  angstroms),  a ,  b7  c, 
a,  /?,  and  y  are  the  cell-edge  lengths  and  angles  of  the  unit  cell, 
respectively. 

We  have  previously  shown1*2*24  that,  when  an  accurate 
intermolecular  potential  is  used,  the  removal  of  the  symmetry 
constraints  in  MP  calculations  has  only  a  very  small  effect  on 
the  final  lattice  energies  and  crystallographic  parameters. 
Moreover,  the  crystal  symmetry,  analyzed  at  the  beginning  and 
at  the  end  of  the  energy  minimization,  remained  unchanged.  In 
the  present  work,  we  have  tested  the  effect  of  removing  the 
imposed  symmetry  constraints  in  a  second  series  of  MP 
calculations  for  17  of  the  51  crystals.  These  calculations  have 
been  done  using  the  algorithm  proposed  by  Gibson  and 
Scheraga25  for  efficient  minimization  of  the  energy  of  a  fully 
variable  lattice  composed  of  rigid  molecules  and  implemented 
in  the  program  LMIN.26  In  these  calculations,  the  parameters  P 
and  Qy  which  specify  the  start  and  the  end  of  the  cubic  feather 
(see  refs  1  and  25  for  details),  were  set  to  20.5  and  20.0, 
respectively. 

As  in  our  earlier  studies,  we  analyzed  the  crystal  symmetry 
of  each  of  the  17  systems  at  the  beginning  and  the  end  of  the 
energy  minimization  to  determine  if  the  crystal  space  group  was 
conserved  in  the  energy  minimization.  The  space  group  is 
considered  to  be  conserved  if  the  symmetry  operations,  as 
defined  in  die  International  Tables  of  Crystallography,27  between 
the  molecule(s)  in  the  asymmetric  unit  cell  and  the  remaining 
molecuie(s)  in  the  unit  cell  are  unchanged  and  if  the  lattice 
parameters  fixed  by  the  lattice  symmetry  have  not  been  modified 
significantly. 

Another  quantity  of  interest  is  the  lattice  energy  of  crystals. 
When  different  crystallographic  phases  exist,  the  lattice  energy 
can  provide  information  about  their  relative  stabilities.  Moreover, 
the  calculated  static  lattice  energy  of  the  crystals  can  be 
compared  to  the  experimental  sublimation  enthalpy  based  on 
the  relation28  —  Atfsubi  =  E  4-  Ko  4-  2 RTy  where  E  is  the  lattice 
energy  and  Ko  is  the  zero-point  energy.  Often,  a  rough  estimate 
of  the  lattice  energy  is  obtained  by  neglecting  the  Ko  term. 
Kitaigorodski18  has  pointed  out  that  considering  the  inaccuracy 
involved  in  the  experimental  determination  of  Atfsubi  and  due 
to  neglect  of  zero-point  energy,  discrepancies  up  to  3—4  kcal/ 
mol  between  the  calculated  and  the  observed  enthalpies  of 
sublimation  are  expected. 

Isothermal— Isobaric  Molecular  Dynamics  Calculations. 
A  more  stringent  test  of  the  intermolecular  interaction  potential 
is  accomplished  through  prediction  of  the  structural  lattice 
parameters  by  isothermal- isobaric  MD  simulations  at  different 
temperatures.  We  have  performed  such  calculations  for  two  of 
the  most  important  energetic  crystals  in  our  database,  namely 
the  monoclinic  phase  of  2,4,6-trinitrotoluene  (TNT)  (entry  14, 
Table  IS  of  the  Supporting  Information)  and  the  tetragonal  phase 


992  J.  Phys.  Cherru  A ,  Vol.  103 ,  Afo.  S,  7999 


Sorescu  et  al. 


f 

(i)  NTR0MA13  (2)  HEVRUV  (3)  JUTGEX 


CH~  NOj 


-CK,  OjN- 


>  NO; 


o^v^A 

9^  NOj 

(6)CUVXOG 


,Vf\" 


q*  Nq 

(7)  CEYDUF 


NO, 


Nq  CH3  Nq  NOj 

(4)  BEOEY  (5)  0Q0BRDC2 

NO,  NOj  Nq 

V  M  ^ 

NO,  Nq 


(1)  FOHMUX  (9)  JOHBUDOl 
CK, 

NqqN 


(ll)VUCBAW  (12)  BECJUO  (13)7NBENZI0  <]4)  TNT(J).  05)  IKTOD 
CHj 

Q,N  Nq  H,C 


s3AvCH’ 

3s@@ 


q*  sqNq 
(IS)  20RHUJ 


qN- 


fa 

0 

H-A 

'  ! 

7  q^ 

wq 

(22)  HASHEO  (23)  KACXEXJ 

no. 


qN  Nq  qN  *q 

(24)  JUWfTW  (25)  UNKUL  (26)  LINJAT 

no.  Na 


C22)  CAXNIY 


wq 

(2X)  BUYPUGJO 


qN  Nq 

(30)  VDCYUJ 


(31)  DATWA1 


Figure  1.  Illustration  of  the  molecules  whose  crystal  structures  were 
studied.  Where  available,  the  corresponding  refcode  entry  in  the 
Cambridge  Structural  Database33  is  indicated. 

of  pentaerythritol  tetranitrate  (PETN  I)  (entry  50,  Table  IS). 
Simulations  were  performed  at  atmospheric  pressure  and 
temperatures  from  100  to  350  K  for  TNT  and  100—400  K  for 
PETN  using  the  algorithm  proposed  by  Nose  and  Klein29  as 
implemented  in  the  program  MDCSPC4B.30  Since  details  of 
the  calculations  are  given  in  ref  1,  we  will  give  only  a  brief 
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Figure  2.  Calculated  percentage  errors  between  the  predicted  and  the 
experimental  values  of  the  lattice  dimensions  a  in  (a),  b  in  (b)  and  c  in 
(c)  for  all  crystals  given  in  Table  IS  of  the  Supporting  Information. 
The  crystal  index  corresponds  to  the  order  number  in  Table  IS. 


Crystal  Index 

Figure  3.  Calculated  structural  shift  factor  F  (eq  3)  for  the  crystal 
structures  in  the  database  as  a  function  of  the  electrostatic  set  of  charges. 
The  crystal  index  corresponds  to  the  order  number  in  Table  1 S  of  the 
Supporting  Information.  The  horizontal  lines  at  1%  and  2%  are  marked 
for  a  more  clear  view  of  the  distribution  of  points. 

description  of  the  computational  parameters.  The  MD  simulation 
cells  consist  of  boxes  containing  16  (2  x  4  x  2)  and  36  (3  x 
3x4)  unit  cells  for  TNT  and  PETN  t  respectively.  The  lattice 
sums  were  calculated  subject  to  minimum-image  periodic 
boundary  conditions  in  all  dimensions.31  The  interactions  were 
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TABLE  1:  Comparison  of  the  Experimental  and  Calculated  Lattice  Energies  for  Different  Sets  of  Electrostatic  Charges  Using 
Molecular  Packing  with  Symmetry  Constraints  8  s 


NTROMA13 

HEVRUV 

JUTGEK 

BECJEY 

QQQBRD02 

CUVXOG 

CEYDUF 

FOHMUK 

JOHBUDOl 

BECJIC 

VUCBAW 

BECJUO 

TNBENZ10 

TNT-phase  I 

TNT-phase  II 

GIMBOT 

JUVNAP 

20RHUJ 

DNTNAP01 

CEDZUG 

HASHAK 

HASHEO 

NACXEU 

JUVMIW 

LINHUL 

UNJAT 

CAXNIY 

BUYPUG10 

TAFDU2 

VDCYUJ 

DAFWAI 

LEIKOA 

DUYREU 

DUYRIY 

DAFGAS 

JAJBEB 

MNPHOL02 

DNOPHLOl 

PICRAC11 

JUPRIV 

ZUGPOG 

ZUGPOJ 

DEFLEF 

PABBOJ 

PABBOJ 

DETDOV 

KUNUM 

VUBZUN 

DEMSOD 

CORYIR 

PERYTN10 

PERYTN01 


—60.10 

-80.84 

-104.16 

-121.86 

-103.18 

-181.68 

-161.05 

-96.24 

-112.26 

-110.73 

-140.86 

-121.83 

-132.55 

-151.98 

-149.40 

-257.38 

-164.15 

-177.78 

-139.38 

-141.29 

-161.73 

-190.42 

-206.61 

-157.15 

-138.27 

-161.74 

-122.41 

-188.16 

-155.25 

-169.39 

-151.70 

-178.37 

-162.77 

-173.64 

-168.90 

-200.03 

-101.37 

-120.42 

-152.55 

-148.23 

-119.24 

-140.83 

-214.53 

-113.31 

—226.63 

-129.86 

-176.93 

-114.56 

-143.62 

-147.60 

-221.09 

-178.65 


lattice  energy  (kJ/rool) 
0.9  HF/6-31G** 

-53.50 

-74.67 

-93.40 

-112.41 

-102.11 

-177.88 

-155.43 

-88.38 
-104.03 
-103.17 
-129.99 
-114.52 
-123.69 
-141.49 
-139.28 
-238.13 
-152.72 
-167.66 
—131.59 
-131.55 
-148.04 
-172.25 
-187.59 
-146.64 
-130.74 
-150.76 
-119.31 
-173.58 
-149.45 
-161.57 
-145.12 
-169.02 
-155.44 
-164.69 
-159.12 
-187.13 
-95.34 
-112.81 
— 141.58 
— 140.90 
— 11 1.60 
-130.12 
-195.65 
-106.49 
-212.98 
-127.15 
-162.62 
-104.71 
-135.33 
-136.87 
-202.46 
-166.87 


B3LYP/6-31G** 

-50.83 
-71.34 
-87.42 
-107.53 
-101.77 
-171.51 
-152.74 
-84.49 
-101.07 
-99.47 
-124.19 
-110.01 
-118.11 
—  135.57 
-133.16 
-225.41 
-146.14 
-160.68 
-126.66 
-128.17 
-142.26 
-163.13 
-177.10 
-141.04 
-127.00 
-144.70 
-117.17 
-164.72 
-145.27 
-156.96 
-141.75 
-162.62 
-151.10 
-159.68 
-153.77 
-180.24 
-94.96 
-110.41 
-13530 
-138.97 
-110.41 
-126.87 
-187.30 
-103.22 
-206.47 
-126.94 
-153.67 
-102.95 
-133.15 
—  132.70 
-193.77 
-161.77 


MP2/6-31G** 

-49.38 
-70.20 
-8531 
-106.08 
-10136 
-168.89 
-151.49 
-81.68 
-98.72 
-96.17 
-121.46 
-107.42 
-114.48 
-130.97 
-128.70 
-21636 
-143.04 
-157.57 
-123.25 
-124.99 
-138.04 
-15831 
-171.47 
-137.90 
-123.95 
-141.01 
-115.95 
-160.95 
.  -143.42 
—  154 25 
-139.03 
-159.09 
-148.62 
-156.96 
-151.05 
-176.83 
-93.61 
-10739 
-129.16 
-135.27 
-108.05 
-124.04 
-179.94 
-10131 
-202.42 
-124.59 
-149.16 
-100.50 
-130.60 
-13134 
-17038 
-156.84 


A#**  (kJ/mol)* 


54.8  ±  4.23* 


70.7  ±  1.7** 


1073  ±  0.636c 
118.4  ±43** 

179.93* 


96.4  ±  1.43* 


91.6  ±2.1*^ 
104.6  ±  4.236g 
105.1  ±  1.6* 


151.9±  2.136aJi 


a  The  corresponding  references  are  given  in  ref  36. 

determined  between  the  sites  (atoms)  in  the  simulation  box  and 
the  nearest-image  sites  within  the  cutoff  distance.  Cutoff 
distances  were  set  at  11.0  A  for  both  crystals.  In  the  initial 
simulation  corresponding  to  the  lowest  temperature,  the  position 
and  orientation  of  the  molecules  in  the  unit  cell  were  taken  to 
be  those  for  the  experimental  structure.  The  initial  velocities  of 
the  centers  of  mass  of  the  molecules  were  selected  at  random, 
but  were  modified  to  eliminate  the  translation  and  rotation  of 
the  bulk  MD  cell.  The  trajectories  were  integrated  for  12  000 
time  steps  (1  time  step  =  2  x  10"15  s),  of  which  2000  steps 
were  equilibration.  In  the  equilibration  period,  the  velocities 
were  scaled  after  every  five  steps  such  that  the  internal 
temperature  of  the  crystal  mimics  the  imposed  external  tem¬ 
perature.  Then,  average  properties  were  calculated  over  the 


remaining  integration  steps  in  the  simulation.  In  subsequent  runs, 
performed  at  successively  higher  temperatures,  die  initial 
configurations  of  the  molecular  positions  and  velocities  were 
taken  from  the  previous  simulation  at  the  end  of  the  production 
run.  The  velocities  were  again  scaled  over  an  equilibration 
period  of  2000  steps,  to  achieve  the  desired  external  temperature, 
followed  by  a  10000-step  production  run. 

Several  types  of  quantities  were  determined  to  obtain 
information  about  structural  parameters  of  the  crystal.  These 
include  the  mean  lattice  geometrical  parameters,  the  cumulative 
mass-center  radial  distribution  functions  (RDF),  and  the  average 
positions  and  orientations  of  the  molecules.  These  quantities 
were  obtained  from  values  calculated  at  every  10th  step  during 
the  trajectory  integrations. 
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TABLE  2:  Lattice  Parameters  and  Energies  Obtained  in  Crystal-Packing  Calculations  without  Symmetry  Constraints8 _ 


final  lattice  parameters 


crystal 

a  (A) 

b  (A) 

c(A) 

a (deg) 

/3  (deg) 

y  (deg) 

final  energy 

HEVRUV 

10.3411  (—0.1) 

103412  (-0.1) 

10.3422  (-0.1) 

89.99  (0.0) 

89.99  (0.0) 

89.99  (0.0) 

-70.22 

BECJEY 

6.2379  (-1.0) 

6.2898  (0.2) 

11.6008  (-2.1) 

100.02  (-0.7) 

80.26  (-0.9) 

118.21  (-0.4) 

-106.11 

TNBENZ10 

9.3251  (-4.6) 

27.2363  (1.1) 

12.8791  (0.4) 

89.99  (0.0) 

89.99  (0.0) 

89.98  (0.0) 

-114.46 

TNT  I 

20.5078  (-3.6) 

6.1176(0.4) 

14.8125  (-1.4) 

90.00  (0.0) 

110.39  (0.2) 

89.99  (0.0) 

-131.02 

TNTH 

14.7733  (-1.5) 

19.2835  (-3.7) 

6.1343  (0.6) 

89.99  (0.0) 

89.99  (0.0) 

90.00  (0.0) 

-128.97 

GIMBOT 

22.3504(0.1) 

5.5723  (0.0) 

14.8260  (1.0) 

90.00(0.0) 

111.13(0.9) 

89.99  (0.0) 

-216.47 

NACXEU 

6.7829  (2.2) 

23.2037  (-0.3) 

7.8659  (0.1) 

90.00  (0.0) 

114.89  (1.5) 

89.99  (0.0) 

-171.46 

JUVMIW 

8.4259  (-1.8) 

10.4800  (-2.0) 

113478  (-1.2) 

90.00  (0.0) 

90.00(0.0) 

89.99  (0.0) 

-137.94 

LINHUL 

6.4673  (-3.7) 

11.9504  (-2.6) 

12.9740  (-2.2) 

63.90  (-03) 

81.62  (-0.2) 

88.69  (-0.6) 

-123.98 

BUYPUG10 

7.7691  (-1.3) 

7.7692  (-1.3) 

10.4845  (-0.64) 

89.99  (0.0) 

90.00(0.0) 

89.99  (0.0) 

-161.00 

DUYRJY 

9.5770(0.1) 

11.5415  (—1.4) 

11.7313  (-1.0) 

76.60  (-0.2) 

89.66  (03) 

88.27  (0.1) 

-156.89 

PICRAC1 1 

9.1385  (-1.3) 

19.0092  (-0.6) 

9.5956  (-1.2) 

90.00  (0.0) 

90.00  (0.0) 

89.99  (0.0) 

-257.50 

JUPRTV 

8.2936  (-0.9) 

16.3874  (-0.7) 

8.6169  (-0.4) 

90.01  (0.0) 

199.58  (1.5) 

89.99  (0.0) 

-133.18 

DEFLEF 

6.8144  (-0.6) 

9.5411  (-1.0) 

18.6262  (0.4) 

90.00(0.0) 

90.00  (0.0) 

89.99  (0.0) 

-180.03 

PABBOJ 

6.6714  (-0.4) 

20.2231  (-0.5) 

11.4702  (-3.3) 

89.99  (0.0) 

90.00(0.0) 

90.00(0.0) 

-201.87 

PETNI 

9.2832  (-1.0) 

9.2841  (-1.0) 

6.5995  (-1.6) 

90.00  (0.0) 

90.00  (0.0) 

90.00  (0.0) 

-170.51 

PETNU 

12.9845  (-2.3) 

13.4333  (-0.4) 

6.8734  (0.6) 

90.00  (0.0) 

90.00(0.0) 

90.00  (0.0) 

-156.88 

a  The  values  in  parentheses  are  the  percent  difference  relative  to  experimental  values. 


TABLE  3:  Predicted  Lattice  and  Volume  Parameters  as  Functions  of  Temperature.  The  Calculated  Thermal  Expansion 
Coefficients  (%)  at  300  K  Are  Also  Indicated. _ _ _ 


T( K) 

a(k) 

b(k) 

c(A) 

a  (deg) 

(deg) 

y  (deg) 

volume  (A3) 

TNT  I 

exptF 

213750 

6.09301 

5.0250 

90.0 

110.14 

90.0 

1828.57 

100.0 

20.5812 

6.13851 

4.8520 

90.00 

11037 

89.99 

1758.71 

200.0 

20.6753 

6.16141 

4.8971 

89.98 

11033 

90.02 

1779.16 

273.1 

20.7553 

*6.18001 

4.9299 

90.01 

110.34 

89.99 

1795.19 

300.0 

20.7839 

6.18751 

4.9452 

89.97 

11031 

90.02 

1802.03 

350.0 

20.8418 

6.20111 

4.9742 

89.99 

11031 

90.00 

181450 

** 

54.0  x  10-*6 

43.4  x  10"6 

34.6  x  10-* 

135.8  x  I0“* 

PETNI 

expd 

9.3776 

93776 

6.7075 

90.0 

90.0 

90.0 

589.85 

100.0 

9.3039 

9.3094 

6.6164 

90.00 

89.99 

89.99 

572.97 

200.0 

93218 

9.3254 

6.6327 

89.98 

89.86 

89.99 

576.48 

273.1 

93399 

93387 

6.6491 

89.98 

89.98 

89.99 

579.91 

300.0 

9.3472 

9.3470 

6.6571 

89.99 

89.98 

90.00 

58150 

350.0 

93576 

93550 

6.6686 

90.01 

90.00 

89.99 

58356 

400.0 

93680 

93720 

6.6756 

89.99 

90.00 

89.99 

585.98 

X 

22.7  x  10-* 

21.2  x  10“* 

31.0  x  10”6 

75.7  x  10~6 

a  The  experimental  values  at  300  K.  *The  units  for  the  liniar  and  volume  expansion  coefficients  are  K”1. 


IV.  Results  and  Discussions  phases  of  the  TNT  and  PETN  crystals.  The  specific  references 

for  all  51  crystals  are  given  in  ref  35. 

Molecular  Packing  Calculations  'with  Symmetry  Con-  The  results  of  MP  calculations  using  the  PCK91  program 

strain ts.  The  51  nitro  compounds  considered  in  this  study  are  ape  presented  in  supplemental  Table  IS.  The  results  in  this  table 

shown  in  Figure  1.  This  set  of  crystals  includes  nitroalkanes,  and  Figure  2  show  that  die  predicted  structural  lattice  parameters 

nitroaromatics,  nitrocubanes,  polymtroadamantanes,  polym-  for  almost  all  of  the  crystals  differ  by  less  than  3%  from  the 

tropolycycloundecanes,  polynitropolycyclododecanes,  hydroxy-  experimental  structures.  The  largest  differences  are  for  die 

nitro  derivatives,  nitrobenzonitriles,  nitrobenzotri azoles,  and  TNBENZ10  (entry  13,  Table  IS  of  the  Supporting  Information) 

nitrate  esters.  Our  selection  includes  some  important  examples  crystal  with  a  maximum  value  of  —4.6%  for  one  of  the  lattice 

of  energetic  crystals  such  as  TNT,  PETN,  and  polynitro  cage  dimensions.  However,  this  decreases  to  —3.6%  when  the  HF 

compounds.  The  structures  of  most  of  these  crystals  have  been  set  of  charges  is  used.  In  addition,  for  the  majority  of  the 

determined  by  X-ray  diffraction  techniques.  Despite  the  gener-  crystals,  there  are  small  rotations  and  practically  no  translations 

ally  poorer  resolution  of  hydrogen  atom  positions  obtained  by  for  the  molecules  in  the  asymmetric  unit  cell.  The  overall 

these  techniques,  we  have  not  done  any  additional  adjustments  accuracy  of  the  predicted  models  is  evident  in  Figure  3,  which 

of  these  positions  to  give,  for  example,  the  standard  bond  shows  the  structural  drift  factors  described  in  eq  3.  For  the 

lengths.32  The  crystal  structures  in  Figure  1  are  denoted  by  the  majority  (69%)  of  crystals,  the  structural  shift  factor  is  between  . 

corresponding  crystal  “refcode”  used  in  the  Cambridge  Struc-  1  and  2,  while  in  23%  of  the  cases  this  factor  is  less  than  1.0. 
tural  Database.33  The  corresponding  names  of  the  molecules  Table  IS  of  the  Supporting  Information  and  Figure  3  show 

are  given  in  ref  34  and  assigned  a  crystal  index  number,  which  the  influence  of  the  set  of  electrostatic  charges  calculated  at 
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refcode.  In  addition,  we  have  studied  different  crystallographic  the  shift  factors  with  the  set  of  electrostatic  charges  employed 
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Crystal  Index 

Figure  4.  Calculated  lattice  energies  in  percentage  difference  relative 
to  the  corresponding  MP2  values.  The  crystal  index  corresponds  to 
the  order  number  in  Table  IS  of  the  Supporting  Information.  The  three 
horizontal  lines  indicate  the  average  deviations  for  the  energies 
calculated  using  the  B3LYP  «pl)  =  2.6%),  0.9*HF  «p2>  =  6.2%) 
and  HF  ((p3)  =  13.6%)  sets  of  charges. 

However,  for  most  crystals,  the  use  of  B3LYP  and  MP2  charges 
improves  the  agreement  between  the  predicted  and  the  experi¬ 
mental  lattice  parameters.  Some  large  improvements  can  be  seen, 
for  example,  in  the  cases  of  JUTGEK  and  VUCBAW  (entries 
3  and  11  of  Table  IS,  respectively),  where  the  relative  errors 
from  the  experimental  lattice  parameters  decrease  by  more  than 
1%  for  some  lattice  dimensions.  However,  in  almost  all  cases, 
the  variations  of  the  geometrical  lattice  dimensions  are  less  than 
1%  with  respect  to  the  change  of  the  set  of  electrostatic  charges. 

We  can  also  see  from  the  results  in  Table  IS  and  Figure  3 
that  when  the  electrostatic  charges  calculated  at  the  HF  level 
are  scaled  by  0.9,  the  predicted  geometrical  parameters  are  very 
close  to  those  obtained  at  the  MP2  level.  Moreover,  the  structural 
shift  factors  appear  generally  to  have  values  intermediate 
between  the  MP2  and  HF  values. 

The  lattice  energies  predicted  by  different  models  are  given 
in  Table  IS  of  the  Supporting  Information.  As  can  be  seen  by 
comparing  the  data  for  MP2,  B3LYP,  and  HF  methods,  the  use 
of  the  correlated  methods  determines  the  decrease  of  the  absolute 
lattice  energy.  This  effect  can  be  understood  as  a  consequence 
of  the  decrease  in  the  absolute  value  of  the  electrostatic 
interaction,  which  has  a  predominant  attractive  character.  The 
variations  in  the  absolute  values  of  the  HF  lattice  energies  are 
between  1.6%  and  30%  relative  to  the  lattice  energies  deter¬ 
mined  at  MP2  level,  with  an  average  difference  of  13.6%  (see 
Figure  4).  The  use  of  the  0.9  scaling  factor  reduces  these 
differences  in  the  range  0.5“  19%  with  an  average  difference 
of  6.2%.  Finally,  the  B3LYP  lattice  energies  are,  as  expected, 
much  closer  to  the  MP2  energies,  with  a  range  of  variation 
between  0.2  and  13.7%  and  with  an  average  difference  of  2.6%. 
These  results  indicate  that  the  lattice  energies  differ  significantly 
for  sets  of  electrostatic  charges  calculated  with  ab  initio  methods 
that  do  or  do  not  include  the  electron  correlation.  These 
differences  can  be  decreased  by  a  factor  of  ~2  when  the  HF 
charges  are  scaled.  Another  important  result  is  that  DFT  methods 
can  provide  charges  that  give  a  similar  accuracy  (within  2.6%) 
for  the  lattice  energy  to  that  determined  at  the  MP2  level.  This 
finding  is  notable  since  the  computational  time  necessary  for 
B3LYP  is  significantly  lower  than  that  for  MP2.  The  corre- 


J.  Phys .  Chem.  A,  Vol.  103 ,  No.  8,  1999  995 

sponding  average  differences  we  have  found  in  the  case  of 
nitramine  crystals4  were  equal  to  12.8%,  4.1%,  and  2.6%, 
respectively,  for  the  HF,  0.9HF,  and  B3LYP  set  of  charges. 
The  coincidence  of  the  results  found  at  the  B3LYP  level  is  an 
indication  that,  in  both  the  previous  and  present  work,  the 
B3LYP  charges  represent  a  good  approximation  and  a  viable 
alternative  to  the  more  computer  time  demanding  MP2  charges. 

In  Table  1  we  compare  the  calculated  lattice  energies  to  the 
available  experimental  sublimation  enthalpies.  Despite  the 
limited  number  of  experimental  values  given  in  Table  1,  it  can 
be  seen  that  a  significant  improvement  in  the  accuracy  of  the 
predicted  lattice  energies  can  be  obtained  by  using  the  electro¬ 
static  charges  determined  by  methods  that  treat  electron  cor¬ 
relation.  The  scaling  of  the  HF  charges  also  leads  to  improve¬ 
ments  in  the  predicted  energies,  but  the  differences  from  the 
experimental  values  are  larger  than  those  obtained  when  the 
charges  are  calculated  with  electron  correlation  methods.  In 
absolute  values,  the  MP2  energies  for  the  majority  of  crystals 
are  within  the  acceptance  range  of  12-17  kJ/mol  of  the 
experimental  lattice  energies,  as  previously  recommended  by 
Kitaigorodsky.18  Exceptions  are  the  QQBRD02  and  GIMBOT 
crystals  where  the  differences  are  larger.  However,  the  lack  of 
agreement  between  the  calculated  and  experimental  lattice 
energies  does  not  correlate  with  the  accuracy  of  the  predicted 
geometrical  parameters.  In  particular,  for  these  two  crystals,  the 
predicted  lattice  parameters  are  quite  good  (see  entries  5  and 
16  in  Table  IS)  with  maximum  structural  shift  factors  of  0.675 
and  0.955,  respectively. 

We  have  also  determined  the  relative  stability  of  some  crystals 
that  have  different  phases.  Specifically,  we  focused  on  the 
relative  stabilities  of  the  polymorphic  phases  of  TNT  and  PETN 
crystals.  In  the  TNT  case,  the  calculated  MP2  lattice  energies 
for  the  monoclinic  and  orthorhombic  phases  are  —130.97  and 
—128.70  kJ/mol,  respectively.  These  values  indicate  that  the 
monoclinic  phase  is  more  stable  than  the  orthorhombic  phase, 
in  agreement  with  experimental  findings.37  Moreover,  the 
difference  between  the  predicted  lattice  energies  of  these  two 
polymorphs,  23  kJ/mol,  represents  the  energy  of  transformation 
from  the  monociinic  to  the  orthorhombic  phases.  This  result 
compares  well  with  the  experimental  value  AH*  =  1.13  kJ/ 
mol.38  In  the  case  of  PETN  crystal,  our  intennolecular  potential 
predicts  that  the  tetragonal  phase  (PERYTN10)  is  more  stable 
than  the  orthorhombic  phase  (PERYTN01),  also  in  agreement 
with  the  experimentally  determined  stability  ranking.39 

Molecular  Packing  Calculations  without  Symmetry  Con¬ 
straints.  The  results  of  molecular  packing  calculations  for  the 
set  of  17  crystals  arbitrarily  chosen  from  the  entire  set  are  given 
in  Table  2.  These  calculations  were  done  using  the  MP2  charges 
only.  As  can  be  seen  by  comparing  the  data  in  Table  2  with 
that  in  Table  IS  of  the  Supporting  Information,  there  is  very 
good  agreement  between  the  geometric  and  energetic  values 
predicted  in  molecular  packing  with  and  without  symmetry 
constraints.  Small  differences  in  the  total  lattice  energies  (<0.5 
kJ/mol)  between  the  constrained  and  unconstrained  calculations 
are  due  to  differences  in  the  evaluation  of  the  dispersion  lattice 
sums.  The  unconstrained  simulations  do  not  use  the  accelerated 
convergence  method  for  this  evaluation.  In  addition,  we  have 
verified  that  the  symmetry  operations  at  the  beginning  and  at 
the  end  of  energy  minimization  are  unchanged.  This  indicates 
that  the  interaction  potential  sufficiently  describes  the  known 
crystallographic  symmetries  of  these  crystals. 

NPT  Molecular  Dynamics  Calculations.  NPT-MD  calcula¬ 
tions  have  been  performed  for  the  most  stable  crystal  phases  of 
TNT  (monoclinic)  and  PETN  (tetragonal,  denoted  PETN  I). 
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Figure  5.  Comparison  of  the  average  fractional  coordinates  and  Euler  angles  of  the  eight  molecules  in  the  unit  cell  of  TNT  (roonoclinic  phase) 
with  the  corresponding  experimental  values. 


These  calculations  used  the  set  of  MP2  charges  only.  The  crystal 
structure  information  resulting  from  NPT-MD  simulations  at 
atmospheric  pressure  and  different  temperatures  is  given  in 
Table  3.  In  both  TNT  and  PETN  I,  the  lattice  dimensions 
obtained  at  T  =  100  K  (Table  3)  are  in  very  good  agreement 
with  those  determined  from  the  MP  calculations  with  symmetry 
constraints  (Table  IS  of  the  Supporting  Information).  This  is 
expected,  since  the  thermal  effects  at  100  K  should  be  minimal 
and  the  thermal  averages  at  this  temperature  should  be  close  to 
the  values  corresponding  to  die  potential  energy  minimum.  At 
300  K,  the  average  lattice  dimensions  of  these  crystals  agree 
very  well  with  the  experimental  values;  the  corresponding 
differences  for  the  cl,  b7  and  c  lattice  dimensions  are -2.30%, 
1.55%,  and  0.53%  for  TNT  and  0.32%,  0.32%,  and  0.75%  for 
PETN  I,  respectively.  In  addition,  in  both  cases  the  angles  of 
the  unit  cell  are  close  to  the  experimental  values.  At  300  K,  the 
difference  between  the  calculated  and  the  experimental  volume 
of  the  unit  cell  is  1.45%  for  TNT  and  1.41%  for  PETN  L 

Figure  5  provides  a  visual  comparison  of  the  average  mass- 
center  ffactionals  and  Euler  angles  for  each  of  the  eight 
molecules  within  die  unit  cell  of  TNT  with  experimental  values. 
Increasing  the  temperature  from  100  to  300  K  does  not  produce 
any  significant  displacement  of  the  molecular  mass-centers  or 
increase  the  degree  of  rotational  disorder.  Similar  conclusions 
about  the  degrees  of  translational  and  rotational  disorder  were 
obtained  in  the  case  of  PETN  I  (not  shown). 

Additional  evidence  for  the  small  degree  of  translation  of 
the  molecules  inside  the  unit  cell  with  increasing  temperature 
can  be  obtained  from  the  mass-center— mass-center  radial 
distribution  functions  (RDFs).  These  are  given  in  Figure  6.  The 
RDFs  for  both  crystals  exhibit  well-ordered  structure,  with 
correlations  at  long  distances  even  at  the  higher  temperatures. 
The  positions  of  the  major  peaks  do  not  change  significantly, 
and  the  main  temperature  effects  are  the  broadening  of  the  peaks 
and  the  partial  overlapping  of  some  of  them.  For  example,  in 


r(A) 


Figure  6.  Variation  of  the  center  of  mass— center  of  mass  radial 
distribution  function  as  a  function  of  temperature  for  TNT  (a)  and  PETN 
<b). 
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the  case  of  PETN  I,  the  centers  of  mass  of  the  two  molecules 
in  the  unit  cell  occupy  the  special  fractional  positions  (0,0,0) 
and  (l/2,l/2,l/2).  There  are  only  three  specific  distances  between 
these  positions  and  they  correspond  to  very  distinct  peaks  in 
the  RDF  spectrum  (see  Figure  6b). 

We  have  also  determined  the  linear  and  volume  expansion 
coefficients  at  300  K  for  these  crystals  using  the  temperature 
dependence  of  the  lattice  dimensions.  The  values  for  these 
quantities  are  given  in  Table  3.  For  TNT,  the  values  of  5.0  x 
10”5  K~!  and  18.0  x  10“5  K_I  have  been  determined  for  the 
linear  and  volume  expansion  coefficients  in  the  temperature 
range  291—318  K.40  The  measured  linear  expansion  coefficients 
are  in  relatively  good  agreement  with  the  values  predicted  by 
our  potential.  Also,  our  calculated  volume  expansion  coefficient 
is  about  28%  smaller  than  the  reported  experimental  value. 

In  the  case  of  PETN  I,  the  values  of  7.65  x  10“5  K~l  and 
23.2  x  1CT5  K”1  have  been  reported  for  the  linear  and  volume 
expansion  coefficients.41  These  values  axe  a  factor  of  3  larger 
than  our  calculated  expansion  coefficients.  Although  the  interac- 
,  tion  potential  can  reasonably  predict  the  individual  lattice 
dimensions  and  volume  at  room  temperature,  these  differences 
suggest  that  it  does  not  adequately  describe  the  magnitude  of 
the  thermal  expansion.  However,  the  thermal  expansion  coef¬ 
ficients  along  the  a  and  b  axes  are  quite  similar,  in  agreement 
with  the  tetragonal  symmetry  of  this  crystal. 

V.  Conclusions 

We  have  investigated  the  degree  of  transferability  of  a  6-exp 
Buckingham  potential  previously  parametrized  using  experi¬ 
mental  information  for  the  a-RDX  crystal1  to  51  non-nitramine 
crystals,  consisting  of  different  types  of  nitroalkanes,  nitroaro- 
matics,  nitrocubanes,  polynitroadamantanes,  polynitropolycy- 
cloundecanes,  polynitropolycyclododecanes,  hydroxynitro  de¬ 
rivatives,  nitrobenzonitriles,  nitrobenzotriazoles,  and  nitrate 
esters.  The  interraolecular  potential  includes  Coulombic  interac¬ 
tions  between  electrostatic  charges.  These  charges  have  been 
determined  from  fits  to  ab  initio  electrostatic  potentials  calcu¬ 
lated  for  the  individual  molecules  in  the  experimental  configura¬ 
tions.  We  have  considered  four  different  electrostatic  models, 
with  charges  determined  at  HF,  B3LYP,  and  MP2  levels  and 
at  the  HF  level  uniformly  scaled  by  a  factor  of  0.9. 

The  tests  of  this  potential  for  the  entire  collection  of  51 
crystals  have  been  performed  using  MP  calculations  with 
symmetry  constraints.  For  a  smaller  set  of  crystals,  we  have 
verified  that  MP  packing  without  symmetry  constraints  predicts 
essentially  the  same  lattice  dimensions  and  energies.  The 
predicted  crystal  structural  parameters  are  in  good  agreement 
with  the  experimental  values  for  most  of  the  crystals,  with 
differences  generally  less  than  3%.  For  92%  of  the  crystals  in 
the  collection,  the  structural  shift  factor  is  less  than  2.0. 

There  is  only  a  small  influence  (generally  less  than  1%)  on 
the  crystallographic  parameters  by  the  set  of  electrostatic  charges 
used.  However,  the  lattice  energies  of  the  crystals  are  signifi¬ 
cantly  influenced  by  the  electrostatic  model.  In  particular,  the 
best  agreement  with  the  experimental  lattice  energies  has  been 
obtained  for  the  MP2  charges.  The  lattice  energies  calculated 
using  the  B3LYP  charges  overestimate  the  MP2  energies  by 
about  2.6%,  while  the  HF  charges  overestimate  the  MP2 
energies  by  13.6%.  The  procedure  of  uniformly  scaling  the  HF 
charges  by  the  0.9  factor  decreases  the  differences  to  about  6.2%. 

In  the  case  of  the  TNT  and  PETN  I  crystals,  the  intermo- 
lecular  potential  describes  the  correct  order  of  stability  of 
different  phases.  Die  predicted  stabilities  monoclinic  >  ortho- 
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rhombic  in  TNT  and  tetragonal  >  orthorhombic  in  PETN  are 
in  accord  with  experimental  findings.57*39 
Moreover,  the  results  of  NPT-MD  simulations  for  ambient 
conditions  of  temperature  and  pressure  support  the  good 
agreement  of  the  predicted  and  experimental  crystallographic 
values. 

The  success  of  the  present  potential  energy  parameters  in 
describing  different  types  of  crystals  containing  molecules  with 
functional  groups  associated  with  explosives  provides  significant 
incentive  to  further  develop  this  model  through  the  incorporation 
of  the  intramolecular  degrees  of  freedom.  This  will  be  done  in 
future  work. 
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